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We study the spherical collapse of a perfect fluid with an equation of state P — kp hy full general 
relativistic numerical simulations. For < fc < 0.036, it has been known that there exists a gen- 
eral relativistic counterpart of the Larson-Penston self-similar Newtonian solution. The numerical 
simulations strongly suggest that, in the neighborhood of the center, generic collapse converges to 
this solution in an approach to a singularity and that self-similar solutions other than this solution, 
■ including a "critical solution" in the black hole critical behavior, are relevant only when the pa- 

^ ' rameters which parametrize initial data are fine-tuned. This result is supported by a mode analysis 

on the pertinent self-similar solutions. Since a naked singularity forms in the general relativistic 
' Larson-Penston solution for < fc ^ 0.0105, this will be the most serious known counterexample 

04 I against cosmic censorship. It also provides strong evidence for the self-similarity hypothesis in gen- 

eral relativistic gravitational collapse. The direct consequence is that critical phenomena will be 
(N . observed in the collapse of isothermal gas in Newton gravity, and the critical exponent 7 will be 

*\ given by 7 ~ 0.11, though the order parameter cannot be the black hole mass. 

o 



o 
o 

(N 



PACS numbers: 04.20.Dw, 04. 25. Dm, 04.40.Nr 



I. INTRODUCTION 



o 

o 

o . 

There is no characteristic scale in general relativity as well as in Newton gravity. A set of field equations is invariant 
^ ' by scale transformations if we assume appropriate matter fields. It implies the existence of scale-invariant solutions to 
the field equations. Such solutions are called self-similar solutions, which are defined by the existence of a homothetic 
Killing vector field. Although the self-similar solutions are only special solutions of Einstein equations, it often has 
. 5^ ' been supposed that these solutions play an important role in situations where gravity is an essential ingredient in a 
. spherically symmetric system (for example, Carr [Q). Such an assumption can be called the self-similarity hypothesis. 
^ ' A spherically symmetric self-similar system of a perfect fluid has been widely researched. Self-similar solutions in 
Newton gravity have been researched in an effort to obtain simple and reafistic solutions of gravitational collapse . 
In particular, the Larson-Penston solution, which is one of the self-similar solutions, is believed to describe the central 
part of generic spherical collapse of isothermal gas. Recent numerical simulations and results of mode analyses strongly 
support this proposition |^-|^ . In general relativity, a spherically symmetric self-similar system was discussed in various 
situations, such as cosmological voids, gravitational collapse, primordial black holes |pl]|-[l4, and so on. The detailed 
structure of self-similar collapse solutions were analyzed The discovery of the black hole critical behavior shed 
light on the importance of a self-similar solution as a critical solution Several very recent works have been 

done in complete classification of self-similar solutions [p6|-p0[ . 

In the context of the black hole critical behavior, the self-similar critical solution is not an attractor. A renormal- 
ization group approach showed that the critical solution has only one repulsive mode pl| , p2[ . The critical exponent 
which appears in the scaling law of the formed black hole mass is equal to the inverse of the eigenvalue of the repulsive 
mode for a perfect fluid case. 
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In the context of cosmic censorship [^Q, a spherical system of a pressureless fluid (dust) has been extensively 
examined since it can be solved exactly. It has been shown that a naked singularity forms in generic spherical collapse 
of dust from an analytic initial density profile |^.^. It is noted that this solution is not self-similar at all. In the 
presence of pressure, self-similar solutions were investigated by Ori and Piran ||l^,|2^,^. For an adiabatic equation 
of state P — kp{0<k< 0.4), they found a discrete set of self-similar solutions which are analytic both at the 
center and at the sonic point. They discovered the general relativistic counterpart of the Larson-Penston solution for 
< A: 0.036. They observed that a naked singularity forms in this solution for < k ^ 0.0105. They also observed 
that there exist analytic self-similar solutions with naked singularity for every k. Harada ||2g[| showed the generic 
occurrence of naked singularity in spherical collapse of a perfect fluid for < A: < 0.01 by numerical simulations using 
the code based on the Hernandez-Misner formulation without the ansatz of self-similarity. 

The aim of this paper is to examine the validity of the self-similarity hypothesis for a spherical system of a perfect 
fluid and to understand the relation of the self-similarity hypothesis, critical behavior and cosmic censorship. In this 
paper, we adopt the geometrized unit. 



II. BASIC EQUATIONS 
A. Einstein equation 

We adopt the comoving coordinates. The line element in a spherically symmetric spacetime is given by 

ds^ = -e'^(*''-)di2 + e"(*''-)dr2 + R^it, r){de^ + sin^ ed<j)^). (2.1) 

It is noted that the comoving coordinates may be able to cover even inside of the apparent horizon. Wc consider a 
perfect fluid 

pt''' ^{p + P)u^'u'' + Pgt'''. (2.2) 
Then the Einstein equations and the equations of motion for the matter are reduced to the following simple form: 
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where m(t,r) is called the Misner-Sharp mass. We assume the following equation of state: 

P - kp, (2.8) 

where we assume that < k < 1. For a barotropic equation of state, the existence of self-similar solution demands the 
above form. Moreover, this equation of state will be valid for isothermal gas in Newton limit and for relativistically 
high-density polytropes Q . We define the following dimensionless functions : 

(2.9) 
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We also define the following zooming coordinates: 



2 



T = -\n{~t), (2.12) 
z^—^. (2.13) 



It is found that Eqs. (2.5) and (|2.6|) can be integrated as 



4/g 2k 



^ a„{t)z^+>'ri i+f' , (2-14) 
e"" = a„(r)7y-TTr5-4^ (2.15) 

where aa{t) and auj{r) are arbitrary functions. This integrabihty is the advantage of the comoving coordinates. These 
arbitrary functions correspond to the freedom of rescahng the time and radial coordinates as t = t{t) and f = r{t). 
Hereafter, we restrict this gauge freedom by choosing — const and Qui — 1. 

Being transformed into the zooming coordinates, the field equations are reduced to 

M + M' ^jjS^{S + S'), (2.16) 
M + M' ^~kT]S^{S + S'), (2.17) 

^ = 1 + a-i ivz-') ^ z^{S + S'f - v^S\S + S'f, (2.18) 



where the derivatives are abbreviated as 

. d 



dr^ (2-19) 
d 

(2.20) 



d\nz 

For later convenience, we define the quantity y which is one third of the ratio of the "average density" of the region 
interior to r to the local density at r, defined as 



If we consider the regular center, then it is found from Eq. (2.3) that 

y^l (2.22) 

at the regular center z = +0. We can define two velocity functions Vz and Vr. The Vz is the velocity of the z = const 
line relative to the fluid element, which is written as 

Vz^-ze^, (2.23) 
while Vr is the velocity of the R = const line relative to the fluid element, which is written as 



B. Self-similar solutions 



For self-similar solutions, we assume that all dimensionless quantities depend only on z, i.e.. 



V = V{^), (2.25) 

S^S{z), (2.26) 

M = M{z), (2.27) 

ct=ct(z), (2.28) 

w = cj(z). (2.29) 
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The field equations are reduced to the following form: 



(In^)'^-^(l-.), 



(1 - vYv,' - {k + yY + (1 + kYr—s-'' 1 - = 0, 



M 



S 



where Vz and Vr are written as 



and 



Vr 



1-y 

k + y' 



(2.30) 
(2.31) 
(2.32) 

(2.33) 
(2.34) 



respectively. It is noted that, although the apparent form of these equations does not seem to be an autonomous 
system, the original system before we have performed explicit integrations is, of course, an autonomous system. We 
can find that the fluid velocity, with respect to i? = const, vanishes only if Vz = or y = 1. The above set of equations 
together with appropriate boundary conditions is enough to determine the unknown functions M{z), S{z) and t](z). 
However, in addition to the above equations, the following dependent equation is used: 



AyV^^-il + kYS-S~^ 



(2.35) 



Using the fact that y — 1/3 is satisfied at z = +0, from Eqs.( 2.30| )-(2.32), we find the behavior of the solution 
around the regular center z = +0 as 



Ch k 2k r / 2(l + 3fc)\ 

M = -^C^D) — z— 1^1 + O (^zWW j 

^1/3, 



2(l + 3fc) 



/ 2(l + 3fc) \ 

l + Olz 3(1+") j 



T] = 2Dz' 



l + O 



and 



where Ck is a constant determined by k as 



/ 2(l + 3fc) \ 
iz 3(1 + '=) j 

/ 2(l + 3fc) \ 

l + O iz j 
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and the parameter D is defined as 



D=- lim z^^r; = 47rp(t, 0)t^ 

2 z^+o 



Therefore, solutions that have regular centers are parametrized by only one parameter D. Now that we find 

at the regular center, we let e'^ be unity at the regular center by choosing the constant as 

a, = {2D)^. 



(2.36) 
(2.37) 
(2.38) 

(2.39) 
(2.40) 
(2.41) 

(2.42) 
(2.43) 



This gauge fixing gives the physical meaning to the parameter D. Then the behavior of Vz around the regular center 
is written as 
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K - -C, 



-2/3 



{2D) + ;23(l + fc) 



1 



2(l + 3fc) 

O ( z 3(1+ 



(2.44) 



The system of equations has an apparent singularity at a point z ~ Zgp at which the relative velocity oi z — const 
world line with respect to the fluid clement is equal to the sound speed, i.e., 



Such a point is called a sonic point. The regularity requires the following condition at the sonic point: 

4yfc - (1 + kfri-^S-'^ = 0. 



(2.45) 



(2.46) 



Every regular solution must cross the sonic point, satisfying Eq. ( ^.46| ). 

Ori and Piran jl^ discovered the band structure of solutions regular both at the center and at the sonic point. 
In particular, there is only a discrete set of solutions which are analytic at the sonic point. Here, we give special 
attention to such solutions. One of such solutions is the flat Friedmann (FF) solution. There are another types of 
analytic solutions which are called the "black hole" type solutions, in which a massive singularity forms at t = and 
after that the mass of the singularity grows linearly with t, and the "repulsive" type solutions, in which the central 
singularity which forms at t = disappears instantaneously and the cloud begins to expand at i = 0. The solutions 
are characterized by the number of oscillations in the velocity fleld . 

Then we consider the behavior of the analytic similarity solutions. For the FF solution, it is found that 



2 

D = - 



1 



3(l + fc)2 



The solution is written as 



Zl + h , 



2 3(l + fc) ^ 



and 



3 V3(l + A:)^ 

■q = 2Dz^, 
1 



14 = -Cj.'(2D) 3(TTfcTz3TrST, 



(2.47) 

(2.48) 

(2.49) 

(2.50) 
(2.51) 

(2.52) 



For the FF, the big crunch occurs at t = 0, i.e., the singularity occurs at the same time everywhere. 

For k < 0.036, there exists a pure collapse solution. It tends to the Larson-Penston solution in Newton gravity in 
the limit k ^ Q. Hereafter we call this solution the general relativistic Larson-Penston (GRLP) solution. For this 
solution, we flnd that the value of the parameter D is given by 

D « 1.439 (2.53) 

for k 0.01. This coincides with the result of Ori and Piran |]l3| , ^ , ^ . 

We have found another two analytic solutions. These solutions are general relativistic counterparts of Newtonian 
self-similar solutions, Hunter (a) and (b) H. We call these solutions GRHA and GRHB solutions, respectively. These 
similarity solutions are displayed in Figs. ]l|(a)-|l|(g). 

The FF is the only solution which has the big crunch singularity. Unlike the FF, the solution of the black hole and 
repulsive types is regular at /: = 0, except for at r = 0. It implies that the dimensionless physical quantities, such as 
M, S and 77, are finite: i.e.. 



M = Moo, 




(2.54) 


S = 'S'oo, 




(2.55) 




Moo 




V = Voc = 




(2.56) 


y = i, 




(2.57) 
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and 



Vr = Vroo. (2.58) 

ai z — ±00. The solution is black hole type if Vroo > and repulsive type if Vroo < 0. It is found that the GRLP 
and GRHB are black hole type solutions while the GRHA is a repulsive type one. From the above equations, we find 

at t = 0. The velocity function Vz diverges as z ^ cxd as 

Vz^-{2D)-rhrj^^S^^z^. (2.60) 

The number of oscillations of Vr, which coincides with the number of zeroes of (y ~ 1) in the domain < z < 00, is 
0, 0, 1 and 2 for the FF, GRLP, GRHA, and GRHB, respectively. The value of D we have obtained for self-similar 
solutions are summarized in Table | for k = 0.001, 0.008, 0.01, and 0.03. 

III. NUMERICAL EVIDENCE FOR CONVERGENCE TO GRLP 
A. Numerical simulations 

In order to see the generic feature of gravitational collapse, we have num er ical ly simulated the spherical collapse of 



a perfect fluid. We have numerically solved the full Einstein equations (2^)-(2/7) by the standard Misner-Sharp code 
without the self-similarity ansatz. The finite difference equations have been given by the staggered-leapfrog scheme. 
The distribution of grid points has been not homogeneous but concentrated in the neighborhood of the center. The 
total number of the grid points has been 10000. See Harada |^ for details of the numerical code and references are 
therein. 

For simplicity we display the results for time symmetric initial data. It should be noted that we have confirmed 
that the results do not change so much for several models in which initial data is not time symmetric. As a set of 
initial data, we have prepared both homogeneous and inhomogeneous balls of a perfect fluid which are momentarily 
static with vacuum external. For the inhomogeneous models, the density profile has been given by 



P = 



Rs. 



Pc,i 

0, Rs.i < R 



< R< Rs,i, ^2 j^-j 



Since we have confirmed that the results do not depend so much on the detailed form of the density profile, the above 
functional form is considered to represent a typical situation. The models which we have simulated are summarized 
in Table ||, where A4 is the Arnowitt-Deser-Misner (ADM) mass of the ball. 

For k — 0.01, we plot the time evolution of the density profile of models A and B in Figs. |^(a) and|^(b), respectively. 
We can see that model B collapses in a self-similar manner near the center as the collapse proceeds. In particular, 
the density profile around the center tends to p oc R^^ in an approach to the occurrence of singularity, which is 
characteristic to the self-similar solutions as we have seen. In order to see more clearly that the collapse approaches 
the self-similar solution, we plot in Figs, ^(a) and ||(b) the time evolution of the density profile of models A and B for 
k = 0.01, respectively. The ordinate and abscissa are dimensionless quantities Awpt'^ and R/{~t), respectively, where 
t is the proper time at the center and t = is chosen as the occurrence of singularity. For comparison, we also plot 
the FF and GRLP in these figures using the relations Airpt^ = {1/2)t]z~^ and R/{—t) = Sz. It is found that model B 



approaches the GRLP while model A approaches the FF in an approach to the singularity. As seen in Figs. |1 

n 



(b) and 

(e), R/{—t) ~ 0.26 and 0.28 at the sonic point for the FF and GRLP, respectively. Therefore, from Figs. 3(a) and 
3(b), it is found that the approach to the FF and GRLP is not only for the subsonic region but also for the supersonic 
region. 

Moreover, in order to see which self-similar solution the collapse approaches, we have calculated the quantity D 
which is defined by 

D = 4np,f, (3.2) 



6 



where pc is the central density. This definition is consistent with Eq. ( ^.41 ). Actually, we have determined the origin 
of t by requiring that the above defined D tends to be constant. The results for k = 0.001, 0.008, 0.01, and 0.03 are 
plotted in Figs. |j(a)-Q(d). Then, we have found that for k = 0.008, 0.01, and 0.03 most models converge to the GRLP. 
Model A for k = 0.001, 0.008, and 0.01 and model C for k — 0.001 turn out to be trivial counterexamples against 
the convergence. We will discuss them later. The results of numerical simulations are summarized in Table The 
resolution of our code has not been sufficient to show the convergence of models B and D to the GRLP for k = 0.001, 
although some tendency towards the GRLP has been observed. For dust collapse (fc = 0), we can easily find that 
the above defined D approaches 2/3, using the Lemaitre-Tolman-Bondi solution. Therefore, from the continuity with 
respect to k, it is expected that the convergence to the GRLP becomes slower as k goes to zero. In fact, since the 
Newtonian approximation becomes good for fc <C 1, it can be said that the convergence to the GRLP for fc <C 1 has 
been confirmed by Newtonian SPH simulations, by Tsuribc and Inutsuka Then, we conclude that the results of 
numerical simulations strongly suggest that generic spherical collapse converges to the GRLP in an approach to the 
singularity occurrence in both space and time. 



B. Interpretation 



As we have seen in Sec. |III A| , most collapse models approach the GRLP, though several models do not approach 
the GRLP. Here we interpret the results analytically. 

First, we consider homogeneous models such as model A for A: = 0.001, 0.008 and 0.01 and model C for k — 0.001. 
For an initially time symmetric homogeneous ball, the evolution of the central region is described by the closed 
Friedmann solution until the rarefaction wave propagates from the surface to the center. The line element of the 
homogeneous central region is written as 

ds^ = -df + a^it) [dx^ + sin^ x{de^ + sin^ ed(j)'^)\ . (3.3) 

The initial value Ui of the scale factor a and the surface value Xs of the comoving coordinate x s-re written using the 
initial density pi and the initial circumferential radius R^j as 

a? = (3.4) 
Xs = Arcsin ^^^^ . (3.5) 

Therefore, the central region of the initially time symmetric ball begins to contract for k > —1/3. We restrict our 
attention to fc > 0. If the sound wave does not reach the center until the central Friedmann region collapses to the 
big crunch singularity, the central region approaches not the GRLP but the FF. For fc <C 1, the central homogeneous 
part is well described by the closed Friedmann solution with dust, which has the parametrized representation as 

l + cos2^ 

a = a, , (3.6) 

t-U^aJe+^sm2ej. (3.7) 



The big crunch occurs at 6* = 7r/2, i.e., t — ti — 7rai/2 = y/Sn/ {32pi). The trajectory of the rarefaction wave with the 
sound speed Cg = Vk which emanates from the surface at t = ti satisfies 

= -Cs, (3.8) 

dt 

which can be integrated as 

X^Xs- 2cs9. (3.9) 
Therefore, the condition for the sound wave not to reach the center before the big crunch is given by 

Xs 
Cs 

Using the free-fall velocity vjf defined as 



> TT. (3.10) 
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we find the condition 



or that the compactness Ai/Rs,i satisfies 




M , 
> -:rk, 



R 



(3.11) 



(3.12) 



(3.13) 



where we have used fc <C 1. Equation ( 3.13 ) can explain the resuhs of the numerical simulations. Condition ( 3.12 ) 
completely agrees with that for Newton gravity derived by Tsuribe and Inutsuka although the present situation 
can be highly relativistic. In particular, the present analysis is valid even for the evolution inside an apparent horizon. 
Although we have discussed the initially time symmetric case, it is easy to derive a similar condition for the collapse 
in which the central homogeneous region which is not swept by the rarefaction wave, is described by the Friedmann 
solution, which may be n ot on ly the closed Friedmann solution but also the flat or open Friedmann solution. The 
result is the same as Eq. (3.13). 

In general, we can enumerate the following trivial counterexamples against the convergence to the GRLP. If the 
central region can be initially described by the Friedmann solution, then the central region does not approach the 
GRLP, but the FF instead, if the big crunch occurs before the rarefaction wave reaches the center. 

Moreover, it is clear that the exact self-similar solutions other than the GRLP do not approach the GRLP. There 
are another kind of counterexamples. We divide the initial data at i? = Rp into two regions: the central region and 
the surrounding region. If the initial data in the central region is the same as those of the exact self-similar solutions 
other than the GRLP, and Rp is so large that the sound wave cannot reach the center until the central singularity 
forms, the collapse in the neighborhood of the center is described not by the GRLP but by the self-similar solution 
initially prepared in the central region. 

Finally, there exists another type of counterexamples, which can be obtained by the exact fine-tuning of parameters 



which characterize the initial data. We will discuss this type of counterexamples in Sec. IV 



Anyway, it is clear that a set of the above trivial counterexamples occupies only zero-measure in the space of the 
whole of regular initial data. 



IV. MODE ANALYSIS 



As we have seen in Sec. [II, the results of numerical simulations suggest that only the GRLP has an attractive 
nature. In order to confirm this, we examine the behavior of modes in linear perturbations of the self-similar solutions. 



A. Perturbation equations 

We consider the spherically symmetric perturbation around the fixed self-similar solution. We attach the suffix o 
for the background solution. Using the rescaling freedom, we set the arbitrary functions and a^^ to the background 
value, i.e., 



We define the perturbation quantities as 



a<T = a^o = (Zf) 1+'= ; 

O-uj = OojO = 1- 



M = Mo(z) [1 + eMi(T, z) + 0{e^)] 
S^Soiz) [l + eSiiT,z) + 0{e^)] , 
7/ = 770(2;) [l + er]i{T,z) + Oie^)] , 
y = yoiz) [1 + eyi{T, z) + 0{e^)] , 



(4.1) 
(4.2) 



(4.3) 
(4.4) 
(4.5) 
(4.6) 
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where e is a small parameter which controls the expansion. Then we find the equations for perturbations up to linear 
order of e as 



M[ = - 



k yi _ 1 
l + k y 1 + k 



(Ml + 



1 



1 



l + k 
^2 



yyi 



^{l + kfyr"^S-\Mt-St) 



l + k 

-4 



y(Mi + fcy-i^i), 



(4.7) 

(4.8) 



(1 - vf 



l + k 



il^ + sA-{l + k){l-y){Sr + S[) 



{k + yf 



1 



l + k 

yi = Ml - r?i - 351, 



m + ^sA + {i + k){k + y)s[ 



(4.9) 
(4.10) 



where we have omitted the suffix o for simpUcity. 

Assuming the time dependence of the perturbative quantities Qx{t,z) = e'^^Qx{z), we find the following set of 
simultaneous equations: 



M[ 



k yi _ A 
l + k y l + k 



{Mi+ky-'Si), 



I 



A 



-yyi 



y{M,+ky-'Si), 



^'^ l + k""' l + k 
(V^ -k){l-y){k + y)y^ 

kV^{l - yf -{k + yf -^^{1 + kfrj-T^yS- 



(1 - 2k) V^{1 - yf - 3k{k + yf + i(l + kfrj-^yS-'' 



+ {l + k)y{V,^l-y) + {k + y)) X] M-_ 
+ 

-{l + k) {V,^{l-y)-k{k + y))X]Sr. 

We can derive another dependent equation as 

{V^^-k){l-y){k + y)y[ 
= [kV,^{l - y) ((1 - y){\nV^)' - 2y{\ny)') - 2{k + y)y(lny)' 



-i(l + kfr]~T-i^yS~^ 



I 



l + k 



(ln77)' + (ln2/)'-4(lnS')' 



+(1 + k)y ((Vf (1 - 2y) + {k + 2y)) (Int/)' + V^{1 - y){\nV^)') A] Mi 

1 , s3 



+ 



fcy/(i - yf -{k + yf - -(1 + kfr^ys- 



+ {l + k)y{V^{l-y) + {k + y))\]M[ 

+ [{I - 2k)V^{l - y) ((1 - y){\nV^)' - 2y{\ny)') - &k{k + y)y{\ny)' 



+i(l + fc)V^?y5-4 



1 - k 



(lnr/)' + (ln2/)'-4(lnS')' 



-(1 + k) {V^{1 - y){\nV^)' - {V^ + k)y{\ny)') A] S, 

+ (1 - 2k) V^{1 - yf - 3k{k + yf + i(l + kfrj'^yS-^ 

-{l + k) (V^{l-y)-k{k + y))X]S[ 

- [V,\l - y){k + y){\nV,^)' + {V^ - k){l -k- 2y)y{\ny)'] y,. 



(4.11) 
(4.12) 



(4.13) 



(4.14) 



where (Iny)' and (InV/)' are given by 
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(Iny)' (InAf)' - (liiry)' - 3{\nSy, (4.15) 
(InK^)' = - 2[^(lnr;)' - 4(ln5)'. (4.16) 

Then, we examine boundary conditions which the perturbation quantities should satisfy at the boundaries. First, 
we consider the regular center z ~ +0. At the regular center, the definition of y implies that the perturbation of y 
must vanish ai z — +0, since the background solution already satisfies the boundary condition in full order. Then we 
obtain 

yi = 0. (4.17) 



at z = +0. From Eq. ( 1.13 ), it implies the following condition: 

Mi+3kSi=0. (4.18) 

Then, the perturbation solutions which are regular at the center for fixed A are parametrized by one parameter A. 
The boundary condition at z = +0 is written as 

2/1 = 0, (4.19) 
^ ^A, (4.20) 

^^-3(TTI)^' ^'-''^ 

where 

A = 77i(0) = ^(t = -l,r = 0). (4.22) 
P 

A only scales yi, Mi, and Si because we are only considering the linear perturbations. Hence, we can set the 
parameter A as A = 1 without loss of generality. 

Next, we consider the sonic point z = z^p. At the sonic point, we require that the dens ity p erturbation is regular. 



It implies that Mi , 5*1 and yi must satisfy the condition that the right-hand side of Eq. ( 4.14 ) vanishes at the sonic 
point. Only for a discrete set of A, there exists a solution of perturbation equations that is regular both at the regular 
center and at the sonic point. Thus we can obtain eigenvalues A and the associated eigenmodes. 

B. Results of mode analysis 

It is found that the system has a gauge mode with the eigenvalue A given by 



The mode functions are given by 



A.i^. (4.23) 



Ml = y(lnM)', (4.24) 

Si^^ilnSy, (4.25) 

yi = y(lnzj)'. (4.26) 

This mode corresponds to the following gauge transformation: 

i-t)^{-t)-e^{-t)^, (4.27) 

r^r (4.28) 

or, equivalently, 



10 



T^T-e—e—\ (4.29) 
A i-k^ 

z^z + e—e~z. (4.30) 



The eigenvalues of physical repulsive modes for k — 0.001, 0.008, 0.01, and 0.03 are summarized in Table [V, where 
A S R is assumed. For the FF and GRLP, there exists no repulsive mode. On the other hand, the GRHA and GRHB 
have one and two repulsive modes, respectively. Therefore, it is found that only the FF and GRLP can describe 
the final stage of the central region of generic collapse. Together with the existence of the "kink" instability in the 
FF and the self-similar solutions which are not analytic at the sonic point ||l^,Q, we conclude that the GRLP is 
the only self-similar solution that can be an attractor. For the GRHA, there exists only one repulsive mode. This 
solution corresponds to the critical solution in the black hole critical behavior. Only when one parameter p, which 
parametrizes initial data, is fine-tuned around the critical value p* for the black hole formation, this solution has 
importance as a critical solution. The critical exponent 7, which appears in the scaling law of the formed black 
hole mass Mbh oc (p — p*)"' ■, is given by the inverse of this lowest repulsive mode. For k — 0.01, the eigenvalue we 
have obtained agrees well with Maison |22). Since the GRHA solution has a repulsive mode, it is not relevant for 
the behavior of generic collapse. In particular, the final stage of the collapse can be described by the GRHA if the 
parameter is exactly fine-tuned, i.e., p = p*. The GRHB has two repulsive modes. It is expected from the mode 
analyses in Newton gravity [Q, that the solution with n oscillations has n repulsive modes. In order for the solution 
with n oscillations to be relevant, n parameters must be fine-tuned. If the fine-tuning is not exact, the perturbation 
grows into nonlinear regime. Then, it is expected that the collapse will approach the GRLP or disperse away. 



V. DISCUSSIONS 



First, we discuss the validity of self-similarity hypothesis. The results of the numerical simulations and mode 
analysis strongly suggest that generic spherical collapse of a perfect fluid with small k converges to the GRLP in an 



approach to the singularity. This means that the GRLP is an attractor solution. Moreover, in Sec. IIIB, we have 
discussed several counterexamples against the convergence to the GRLP. It is surprising that these counterexamples 
are exactly self-similar in the neighborhood of the center or at least asymptotically self-similar. It should be noted 
that non-flat Friedmann solution also approaches the FF asymptotically in an approach to the big crunch. Therefore, 
we can conjecture that any cloud of a perfect fluid collapses in a self-similar manner in an approach to the singularity. 

Next we discuss the implications of the convergence to the GRLP in the context of the cosmic censorship. The 
cosmic censorship conjecture states that a naked singularity does not form in the gravitational collapse which develops 
from generic initial data with matter fields which obey a physically reasonable equation of state. For spherical collapse, 
the convergence to the GRLP, which we have observed in this paper, means that a naked singularity forms in generic 
collapse for an equation of state P ~ kp ior < k ^ 0.0105, because the GRLP has a naked singularity for that range 
of fc 11,11,11 . 

Here we should give the precise terminology of a naked singularity. In this article, we refer to a singularity that 
can be seen by some observer as a naked singularity. In contrast, a naked singularity which can be seen from infinity 
is called a globally naked singularity. Whether a naked singularity is globally naked is determined not only by the 
central region but also by the surrounding region. In fact, a naked singularity can be globally naked through the 
matching of the central region with an appropriate surrounding region. 

Now that we have the precise terminology, we can discuss the consistency of our results with previous works on 
the black hole critical behavior. At first sight, our results seem to be inconsistent with the formation of an apparent 
horizon observed in numerical simulations showing the black hole critical behavior. In fact, this is not the case. Since 
the convergence is only for the neighborhood of the center, we cannot say whether the formed naked singularity is 
locally naked or globally naked. Because the formation of an apparent horizon only implies the existence of an event 
horizon outside or coinciding with it, it does not exclude the formation of locally naked singularity at the center. 

If the cosmic censorship is true, then there are three possibilities. One is that deviations from spherical symmetry 
may play a crucial role in the nakedness of the formed singularity. Although there has been no systematic study on the 
effect of violation of spherical symmetry in inhomogeneous gravitational collapse, Shapiro and Tcukolsky [ pl| reported 
some numerical results that suggest the occurrence of naked singularity in the axisymmetric collapse of coUisionless 
particles. In contrast, Iguchi et al. [p2|-p^ and Nakao et al. |^ reported some kind of instability along the Cauchy 
horizon associated with a globally naked singularity. The second possibility is that the small value of k is not allowed 
for extremely high-density matter fields. However, it seems to be strange that the consistency of classical theory of 
gravity restricts the equation of state for high-density matter fields, which is determined by a collection of various 
microscopic physics. The third possibility is that the fluid description for high-density matter may be responsible. 
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Whether or not the cosmic censorship conjecture is true, the convergence to the GRLP strongly suggests that there 
can appear an extremely high-density or high-curvature region which can be seen by an observer. Even for such an 
"approximate" naked singularity, it has been shown that explosive radiation is emitted due to quantum effects [^6|-^3| . 
Furthermore, in a practical sense, if the curvature scale reaches the Planck scale, it should be regarded as a singularity 
because it is considered beyond the scope of classical general relativity. 

Self-similar solutions we have obtained here approach those in Newton gravity in the limit A; — > 0. Therefore, 
the important consequence is that critical phenomena associated with the Hunter (a) Newtonian self-similar solution 
should be observed in the collapse of isothermal gas in Newton gravity. These critical phenomena will be very similar 
to the critical phenomena in the black hole formation in general relativity. Only one parameter p has to be fine-tuned 
closely to the critical value p* . In particular, some order parameter A follows the scaling law ^ cx {p — p*)'^ in the 
near critical regime, where 7 is given by the inverse of the eigenvalue of the only one repulsive mode of the Hunter 
(a) solution. Unfortunately, the eigenvalue of the repulsive mode of the Hunter (a) solution has not been known yet. 
However, by extrapolating our results on the GRHA to the limit fc — > 0, we can predict that the critical exponent 7 
is given by 7 « 0.11. The candidate for the order parameter A will be, for example, the mass of the initially formed 
core, if we assume the realistic equation of state for dense gas. 



VI. CONCLUSIONS 



The results of the numerical simulations and mode analysis strongly suggest that the general relativistic Larson- 
Penston solution is an attractor solution of spherically symmetric gravitational collapse of a perfect fluid with an 
adiabatic equation of state P = kp for Q < k 0.036 in general relativity. Since a naked singularity forms in 
the general relativistic Larson-Penston solution for < fc ^ 0.0105, the analysis in this paper means the violation 
of cosmic censorship in spherically symmetric case. This will be the strongest known counterexample against the 
cosmic censorship ever. This also provides a strong evidence for the self-similarity hypothesis in general relativistic 
gravitational collapse. 
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TABLE I. Values of D for self-similar solutions. 





k = 0.001 


0.008 


0.01 


0.03 


FF 


2/(3 X 1.001^) 


2/(3 X l.OOS'') 


2/(3 X l.Ol'') 


2/(3 X 1.03^) 


GRLP 


1.640 


1.480 


1.439 


1.119 


GRHA 


1.675 X 10-^ 


1..345 X 10-^ 


1.26.5 X 10^ 


7.204 X 10^ 


GRUB 


7.171) X iO* 


1.91):! X iO' 


i. ill X iO* 


l.()5l) X iO* 



TABLE IL Models for numerical simulations. 



Models 


Initial Density Profile 


Initial Compactness {M./Rs,i) 


A 


Homogeneous 


1/10 


B 


Inhomogeneous 


1/10 


C 


Homogeneous 


1/30 


D 


Inhomogeneous 


1/30 



TABLE HI. Asymptotic behavior of the collapse models. 



Models 


k = 0.001 


0.008 


0.01 


0.03 


A 


FF 


FF 


FF 


GRLP 


B 


GRLP? 


GRLP 


GRLP 


GRLP 


C 


FF 


GRLP 


GRLP 


(Dispersion) 


D 


GRLP? 


({RLl' 


(;rli' 


(DisiK-rsitHi) 



TABLE IV. Eigenvalues A for repulsive modes of self-similar solutions. 





A- = O.OOi 


O.ODS 


D.Oi 


().();', 


FF 


None 


None 


None 


None 


GRLP 


None 


None 


None 


None 


GRHA 


9.39 


8.88 


8.75 


7.62 


GRHB 


5.43 


5.10 


5.02 


4.27 




5.62 X 10 


4.90 X 10 


4.71 X 10 


3.30 X 10 



14 




15 




16 



1e+06 



10000 



100 



0.01 



0.0001 



FF 
GRLP 
GRHA 
GRHB 




17 



(g) 
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R/(-t) 



10 1 00 1 000 



FIG. 1. Self-similar solutions for k = 0.01. (a) InM, (b) \nS, (c) Inr], (d) y, (e) -14, and (f) -Vr are plotted. In (g), the 
ordinate and abscissa are 47rpt^ and R/{—t), respectively. 
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0.0001 



FIG. 2. Time evolution of the density profile for models (a) A and (b) B for A; = 0.01 are plotted. The ordinate and abscissa 
are the density p and the circumferential radius R, respectively. The unit is chosen so that the ADM mass M is unity. 
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1e-14 I ' ' ' ' 1 

1e-2 1 1e+2 1e+4 

(b) ^/(-^> 
FIG. 3. Time evolution of the density profile for models (a) A and (b) B for A; = 0.01. The ordinate and abscissa are A-irpt^ 
and Rj{—t), respectively. For comparison, the FF and GRLP are also plotted. 
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'2 for (a) k 
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FIG. 4. D = Anpct^ for (a) k = 0.001, (b) k = 0.008, (c) k = 0.01, and (d) k = 0.03 are plotted. The abscissa is the central 
density pc. The values of D for the FF {D = 0.6653, 0.6561, 0.6535 and 0.6284) and for the GRLP {D = 1.640, 1.480, 1.439 
and 1.119) for k = 0.001, 0.008, 0.01 and 0.03 are also denoted, respectively. 



22 



